rm(list = ls())
library(foreign)
library(fwb)
library(dplyr)
library(brglm2)
library(beepr)
library(modelsummary)
library(lmtest)
library(sandwich)
library(marginaleffects)
library(kableExtra)
library(marginaleffects)
library(ggplot2)

options(modelsummary_output = "latex_tabularray")
# Remove the current version first (optional, but recommended)
remove.packages("modelsummary")
set.seed(555)

# Install the older version from the CRAN archives
remotes::install_version("modelsummary", version = "0.10.0")
shame2 = read.csv("/Users/cbrandt/Dropbox/AAA_Projects/UNSC/R code/replication_data/brandt_kotajoki_naming_and_shaming_at_the_UNSC_panel_data.csv")

universal_r = 1000
shame2$e10 = ifelse(shame2$p5==1, 0, 1 )
shame2$shamed_dummy <- as.integer(shame2$shamed_dummy)

cm  = c("cw"=  "Civil War",  "isfighting3"="Fighting NSA",
        'ext_confirmed_lag' = 'External Support',
        'ext_confirmed_lag:p5' = 'External Support*P5',
        'p5:ext_confirmed_lag' = 'External Support*P5',
        'best_fatality_estimate_lag_log' = 'One-Sided Violence',
        'best_fatality_estimate_lag_log:p5' = 'One-Sided Violence*P5',
        'p5:best_fatality_estimate_lag_log' = 'One-Sided Violence*P5',
        'ext_confirmed_lag:best_fatality_estimate_lag_log' = 'Ext. Sup.*OSV',
        'best_fatality_estimate_lag_log:ext_confirmed_lag' = 'Ext. Sup.*OSV',
        "no_of_allies_fighting_rebel_dummy"= "Allies",
        "no_of_allies_fighting_rebel_dummy:p5"= "Allies*P5",
        "p5:no_of_allies_fighting_rebel_dummy"= "Allies*P5",
        'ext_confirmed_lag:best_fatality_estimate_lag_log:p5' = 'Ext. Sup.*OSV*P5',
        'best_fatality_estimate_lag_log:ext_confirmed_lag:p5' = 'Ext. Sup.*OSV*P5',
        "best_fatality_estimate_lag_log:no_of_allies_fighting_rebel_dummy" = "Allies*OSV",
        'no_of_allies_fighting_rebel_dummy:best_fatality_estimate_lag_log' = 'Allies*OSV',
        
        "ext_nsa_confirmed_lag"="Ext Support to NSAG",
        "ext_nsa_confirmed_lag:p5"="Ext Support to NSAG*P5",
        "p5:ext_nsa_confirmed_lag"="Ext Support to NSAG*P5",
        
        "post911"="post911",
        "any_osv" = "Any OSV",
        "p5" = "P5", 
        "v2x_clphy_lag" = "Physical Integrity",
        "v2x_clphy_lag:best_fatality_estimate" = "One-Sided Vio. x Phys. Int.",
        "any_osv:ext_confirmed" = "Ext Support*Any OSV",
        
        "polity2" = "Regime Type",
        "battledeaths_lag_log"= "Intensity",
        "strength" = "Strength",
        "incompatibility"= "Secessionist",
        "hrw_prev_2" = "SV: HRW2",
        "state_prev_2" = "SV: State2",
        "ai_prev_2" = "SV: AI2",
        "rebpolwing" = "Political Wing",
        "pk_lag_agg_log" = "Peacekeeping",
        
        
        "prior_no_russia" = "Prior Shame*",
        "prior_no_china" = "Prior Shame*",
        "prior_no_us" = "Prior Shame*",
        "prior_no_uk" = "Prior Shame*",
        "prior_no_france" = "Prior Shame*",
        "prior_no_russia_china" = "Prior Shame*",
        "prior_no_western" = "Prior Shame*",
        
        
        "any_member_shamed_year_prior"= "Prior Shame",
        "total_resolutions_log_lag" = "Resolutions",
        "t_a"="Terrorism1",
        "t_b"="Terrorism2",
        "k_a"="Terrorism3",
        "k_b"="Terrorism4",
        
        "us_shamed_same_year_prior"= "US Shamed Prior",
        "russia:us_shamed_same_year_prior" = "US Shamed Prior*Russia Fe",
        "us_shamed_same_year_prior:russia" = "US Shamed Prior*Russia Fe",
        "china:us_shamed_same_year_prior" = "US Shamed Prior*China FE",
        "us_shamed_same_year_prior:china" = "US Shamed Prior*China FE",
        "uk_shamed_same_year_prior"= "UK Shamed Prior",
        "france_shamed_same_year_prior"= "France Shamed Prior",
        "russia_shamed_same_year_prior"= "Russia Shamed Prior",
        "china_shamed_same_year_prior"= "China Shamed Prior",
        
        "western" = "Western",
        "russ_china" = "Rus/China",
        "western_shame_prior" = "Western Shamed Prior", 
        "russ_china_shame_prior" = "Russia/China Shamed Prior",
        
        "western_shame_prior:russ_china" =  "Western Prior*RC", 
        "western_shame_prior:russia" =  "Western Prior*Russia", 
        "western_shame_prior:china" =  "Western Prior*China", 
        "western_shame_prior:western" =  "Western Prior*Western",
        "russ_china_shame_prior:western" = "RC Prior*Western",
        "p5_shamed_same_year_prior"= "P5 Shamed Prior",
        "p5:p5_shamed_same_year_prior" = "P5 Shamed Prior*P5",
        "e10_shamed_same_year_prior"= "E10 Shamed Prior",
        "e10_shamed_same_year_prior:p5" = "E10 Shamed Prior*P5",
        "rolling_shaming_group_sum_5yr" = "Rolling by Group", 
        "rolling_shaming_state_sum_5yr" = "Rolling by State",
        "nego" = "Negotiations",
        "usa" = "US FE",
        "russ_china_shame_prior:usa" = "Russia/China Prior*US FE",
        
        "uk" = "UK FE",
        "russ_china_shame_prior:uk" = "Russia/China Prior*UK FE",
        "france" = "France FE",
        "russ_china_shame_prior:france" = "Russia/China Prior*Fr FE",
        "russia" = "Russia FE",
        "china" = "China FE",
        
        "usa:russ_china_shame_prior" = "Russia/China Prior*US FE",
        "uk:russ_china_shame_prior" = "Russia/China Prior*UK FE",
        "france:russ_china_shame_prior" = "Russia/China Prior*Fr FE",
        
        "us_shamed_same_year_prior:usa" = "US Shamed Prior*US FE",
        "usa:uk_shamed_same_year_prior" = "UK Prior*US FE",
        "usa:france_shamed_same_year_prior" = "France Prior*US FE",
        "usa:russia_shamed_same_year_prior" = "Russia Prior*US FE",
        "usa:china_shamed_same_year_prior" = "China Prior*US FE",
        
        "uk:us_shamed_same_year_prior" = "US Shamed Prior*UK FE",
        "uk:uk_shamed_same_year_prior" = "UK Prior*UK FE",
        "uk:france_shamed_same_year_prior" = "France Prior*UK FE",
        "uk:russia_shamed_same_year_prior" = "Russia Prior*UK FE",
        "uk:china_shamed_same_year_prior" = "China Prior*UK FE",
        
        "france:us_shamed_same_year_prior" = "US Shamed Prior*Fr FE",
        "france:uk_shamed_same_year_prior" = "UK Prior*Fr FE",
        "france:france_shamed_same_year_prior" = "France Prior*Fr FE",
        "france:russia_shamed_same_year_prior" = "Russia Prior*Fr FE",
        "france:china_shamed_same_year_prior" = "China Prior*Fr FE",
        
        
        "russia:uk_shamed_same_year_prior" = "UK Prior*Ru FE",
        "russia:france_shamed_same_year_prior" = "France Prior*Ru FE",
        "russia:russia_shamed_same_year_prior" = "Russia Prior*Ru FE",
        "russia:china_shamed_same_year_prior" = "China Prior*Ru FE",
        
        "china:uk_shamed_same_year_prior" = "UK Prior*Ch FE",
        "china:france_shamed_same_year_prior" = "France Prior*Ch FE",
        "china:russia_shamed_same_year_prior" = "Russia Prior*Ch FE",
        "china:china_shamed_same_year_prior" = "China Prior*Ch FE",
        "ally_shamed_prior_2_years" = "ally Shamed Prior",
        'strength_lag' = "Strength",
        "t_a_lag"="Terrorism1",
        "t_b_lag"="Terrorism2",
        "k_a_lag"="Terrorism3",
        "k_b_lag"="Terrorism4",
        'neg_lag' = 'Negotiations',
        
        'islam' = "Islamist",
        'post_911' = "Post 9/11",
        'islam:post_911' = "Islamist*Post 9/11",
        
        "year" = "Time",
        
        "e10_shamed_same_year_prior:best_fatality_estimate_lag_log" = "E10 Shamed Prior*OSV",
        
        "best_fatality_estimate_lag_log:p5_shamed_same_year_prior" = "P5 Shamed Prior*OSV",
        
        "e10_shamed_same_year_prior:best_fatality_estimate_lag_log:p5" = "E10 Shamed Prior*OSV*P5",
        
        "best_fatality_estimate_lag_log:p5:p5_shamed_same_year_prior" ="P5 Shamed Prior*OSV*P5", 
        
        "e10_shamed_same_year_prior:ext_confirmed_lag" = "E10 Shamed Prior*Ext. Sup.",
        
        "ext_confirmed_lag:p5_shamed_same_year_prior" = "P5 Shamed Prior*Ext. Sup.",
        "e10_shamed_same_year_prior:ext_confirmed_lag:p5" = "E10 Shamed Prior*Ext. Sup.*P5",
        "ext_confirmed_lag:p5:p5_shamed_same_year_prior" = "P5 Shamed Prior*Ext. Sup.*P5",
        "(Intercept)"="Intercept"
)

#######Main Analysis ======
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior

formula2 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5  +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag 

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = ~unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = ~unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = ~unsc_gwn, R = universal_r)



# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c( "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title = "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/main_analysis.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#marginal effects of main analysis ------
# --- 1. Define the Specific Points of Interest ---
# Pre-calculate the "mean + 1 standard deviation" value for clarity
mean_plus_sd <- mean(shame2$best_fatality_estimate_lag_log, na.rm = TRUE) + 
  sd(shame2$best_fatality_estimate_lag_log, na.rm = TRUE)

# Create a vector with the two exact values you want
fatality_points <- c(1, mean_plus_sd)

# --- 2. Calculate Predicted Probabilities at Those Points ---
# The datagrid now uses the 'fatality_points' vector
pred_specific_points <- predictions(
  model3,
  newdata = datagrid(
    p5 = c(0, 1),
    best_fatality_estimate_lag_log = fatality_points
  )
)

# You can print this object to see the four predicted probabilities in a table
print(pred_specific_points)

p_low_fatal_e10 <- pred_specific_points %>% filter(best_fatality_estimate_lag_log == 1&p5==0) %>% pull(estimate)
p_high_fatal_e10    <- pred_specific_points %>% filter(best_fatality_estimate_lag_log == mean_plus_sd&p5==0) %>% pull(estimate)

percent_increase_fatalities_e10 <- ((p_high_fatal_e10 - p_low_fatal_e10) / p_low_fatal_e10) * 100
percent_increase_fatalities_e10


p_low_fatal_p5 <- pred_specific_points %>% filter(best_fatality_estimate_lag_log == 1&p5==1) %>% pull(estimate)
p_high_fatal_p5    <- pred_specific_points %>% filter(best_fatality_estimate_lag_log == mean_plus_sd&p5==1) %>% pull(estimate)

percent_increase_fatalities_p5 <- ((p_high_fatal_p5 - p_low_fatal_p5) / p_low_fatal_p5) * 100
percent_increase_fatalities_p5



# --- 1. Calculate Predicted Probabilities ---
# The datagrid now creates a 2x2 grid of all combinations of p5 and ext_confirmed_lag.
pred_confirmation_interaction <- predictions(
  model3,
  newdata = datagrid(
    p5 = c(0, 1),
    ext_confirmed_lag = c(0, 1)
  )
)

# Print the table of the four predicted probabilities
print(pred_confirmation_interaction)

# Pull the two values
p_no_support <- pred_confirmation_interaction %>% filter(ext_confirmed_lag == 0&p5==1) %>% pull(estimate)
p_support    <- pred_confirmation_interaction %>% filter(ext_confirmed_lag == 1&p5==1) %>% pull(estimate)

# Compute percent increase
percent_increase <- ((p_support - p_no_support) / p_no_support) * 100
percent_increase

